EDS 223: week 1 lab
Ruth Oliver
2022-09-28
introduction
In this lab, we’ll explore the basics of map-making in R using the tmap package. The following materials are modified from Chapter 9 of Geocomputation with R by Rovin Lovelace
prerequisites
install.packages('sf')
install.packages('raster')
install.packages('dplyr')
install.packages('spData')
install.packages('tmap')
install.packages('leaflet')
install.packages('ggplot2')
install.packages('gifski')library(sf)
library(raster)
library(dplyr)
library(spData)
library(spDataLarge)
library(tmap)map-making basics
Let’s start with a pre-loaded spatial object representing the states of New Zealand
nz## Simple feature collection with 16 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 1090144 ymin: 4748537 xmax: 2089533 ymax: 6191874
## Projected CRS: NZGD2000 / New Zealand Transverse Mercator 2000
## First 10 features:
## Name Island Land_area Population Median_income Sex_ratio
## 1 Northland North 12500.561 175500 23400 0.9424532
## 2 Auckland North 4941.573 1657200 29600 0.9442858
## 3 Waikato North 23900.036 460100 27900 0.9520500
## 4 Bay of Plenty North 12071.145 299900 26200 0.9280391
## 5 Gisborne North 8385.827 48500 24400 0.9349734
## 6 Hawke's Bay North 14137.524 164000 26100 0.9238375
## 7 Taranaki North 7254.480 118000 29100 0.9569363
## 8 Manawatu-Wanganui North 22220.608 234500 25000 0.9387734
## 9 Wellington North 8048.553 513900 32700 0.9335524
## 10 West Coast South 23245.456 32400 26900 1.0139072
## geom
## 1 MULTIPOLYGON (((1745493 600...
## 2 MULTIPOLYGON (((1803822 590...
## 3 MULTIPOLYGON (((1860345 585...
## 4 MULTIPOLYGON (((2049387 583...
## 5 MULTIPOLYGON (((2024489 567...
## 6 MULTIPOLYGON (((2024489 567...
## 7 MULTIPOLYGON (((1740438 571...
## 8 MULTIPOLYGON (((1866732 566...
## 9 MULTIPOLYGON (((1881590 548...
## 10 MULTIPOLYGON (((1557042 531...- we’re going to start by using the tmap package to make some basic maps
- tmap can work with spatial objects of a variety of clasess, meaning it’s highly versatile
- this approach relies on a series of functions that typically start with “tm_”
the first element is always “tm_shape”
tm_shape(nz) +
tm_fill()
now let’s plot just the boundaries
tm_shape(nz) +
tm_borders()
and the shapes and boundaries together
tm_shape(nz) +
tm_fill() +
tm_borders()
map objects
- tmap can store maps as objects
- this means that we store a base map and add additional layers later
map_nz <- tm_shape(nz) +
tm_polygons()
class(map_nz)## [1] "tmap"- we can add new shapes on top of the base map
- when we add a new shape, all subsequent aesthetic functions refer to it, until we add a new shape
In this case, we’re adding a layer with information on elevation and this layer to have 70% transparency.
map_nz1 <- map_nz +
tm_shape(nz_elev) +
tm_raster(alpha = 0.7)
map_nz1 ## stars object downsampled to 877 by 1140 cells. See tm_shape manual (argument raster.downsample)
we can add points designating high points in the country
map_nz2 <- map_nz1 +
tm_shape(nz_height) +
tm_dots()
map_nz2## stars object downsampled to 877 by 1140 cells. See tm_shape manual (argument raster.downsample)
aesthetic basics
- up until now, we’ve been working with the default aesthetics
- there are 2 types of aesthetics: fixed and those that change with the value of a variable
- tmap works differently than ggplot2 and doesn’t use the “aes()” function
Let’s start by changing some fixed aesthetics…First, let’s change the color used to fill the NZ shapes.
tm_shape(nz) +
tm_fill(col = "red")
now change the color of the boundaries
tm_shape(nz) +
tm_borders(col = "blue")
and the width of the boundary lines
tm_shape(nz) +
tm_borders(lwd = 3)
and the line type of the boundary lines
tm_shape(nz) +
tm_borders(lty = 2)
all together now!
tm_shape(nz) +
tm_fill(col = "red", alpha = 0.3) +
tm_borders(col = "blue", lwd = 3, lty = 2)
Now let’s change the colors based on a value. We noticed that the New Zealand dataset has a column with each state’s land area
nz## Simple feature collection with 16 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 1090144 ymin: 4748537 xmax: 2089533 ymax: 6191874
## Projected CRS: NZGD2000 / New Zealand Transverse Mercator 2000
## First 10 features:
## Name Island Land_area Population Median_income Sex_ratio
## 1 Northland North 12500.561 175500 23400 0.9424532
## 2 Auckland North 4941.573 1657200 29600 0.9442858
## 3 Waikato North 23900.036 460100 27900 0.9520500
## 4 Bay of Plenty North 12071.145 299900 26200 0.9280391
## 5 Gisborne North 8385.827 48500 24400 0.9349734
## 6 Hawke's Bay North 14137.524 164000 26100 0.9238375
## 7 Taranaki North 7254.480 118000 29100 0.9569363
## 8 Manawatu-Wanganui North 22220.608 234500 25000 0.9387734
## 9 Wellington North 8048.553 513900 32700 0.9335524
## 10 West Coast South 23245.456 32400 26900 1.0139072
## geom
## 1 MULTIPOLYGON (((1745493 600...
## 2 MULTIPOLYGON (((1803822 590...
## 3 MULTIPOLYGON (((1860345 585...
## 4 MULTIPOLYGON (((2049387 583...
## 5 MULTIPOLYGON (((2024489 567...
## 6 MULTIPOLYGON (((2024489 567...
## 7 MULTIPOLYGON (((1740438 571...
## 8 MULTIPOLYGON (((1866732 566...
## 9 MULTIPOLYGON (((1881590 548...
## 10 MULTIPOLYGON (((1557042 531...Let’s try to plot the Land_area column. We might think that the following works, but it doesn’t!
#tm_shape(nz) +
# tm_fill(col = nz$Land_area)Instead, tmap is expecting a character string naming the attribute associated with the geometry
tm_shape(nz) +
tm_fill(col = "Land_area")
We can also add an argument that updates the title of the legend
tm_shape(nz) +
tm_fill(col = "Land_area", title = "Area")
We can even make it more precise using the “expression” function
tm_shape(nz) +
tm_fill(col = "Land_area", title = expression("Area (km"^2*")")) +
tm_borders()
color settings
- how we set and define colors can radically change the interpretation of our maps
- the style argument has several options for breaking data into bins
- style = “pretty” (default) rounds breaks into evenly spaced whole numbers, where possible
- style = “equal” divides input values into bins with equal range (best for uniform distributions)
- style = “quantile” puts the same number of observations into each bin
- style = “jenks” identifies groups of similar values and maximizes differences between bins
Notice how the following maps display the same data, but look quite different
tm_shape(nz) +
tm_polygons(col = "Median_income", style = "pretty")
tm_shape(nz) +
tm_polygons(col = "Median_income", style = "equal")
tm_shape(nz) +
tm_polygons(col = "Median_income", style = "quantile")
tm_shape(nz) +
tm_polygons(col = "Median_income", style = "jenks")
We can also define custom bins
breaks = c(0, 3, 4, 5) * 10000
tm_shape(nz) +
tm_polygons(col = "Median_income", breaks = breaks)
- in some cases we might not want to use bins
- style = “cont” displays colors as a continuous spectrum
- style = “cat” uses a unique vale for each categorical value
map_nz +
tm_shape(nz_elev) +
tm_raster(alpha = 0.7,
style = "cont")## stars object downsampled to 877 by 1140 cells. See tm_shape manual (argument raster.downsample)
map_nz +
tm_shape(nz) +
tm_polygons(col = "Island",
style = "cat")
map layout
- now that we have the basics, we can turn to all the other elements that make a cohesive map
- tmap has lots of options, but we will explore just a few
To clearly give readers the context of our map, we can include a compass and scale bar
map_nz +
tm_compass(type = "4star", position = c("left", "top")) +
tm_scale_bar(breaks = c(0, 100, 200), text.size = 1)
Instead of using a compass and scale bar, we could add latitude/longitudes graticules
map_nz +
tm_graticules()
We can also update the background color
map_nz +
tm_graticules() +
tm_layout(bg.color = "lightblue")
faceted and animated maps
- we might have data that varies over time and we want to look at the how it changes
- one approach is by using faceted plots, a series of plots
urb_1970_2030 <- urban_agglomerations %>%
filter(year %in% c(1970, 1990, 2010, 2030))
tm_shape(world) +
tm_polygons() +
tm_shape(urb_1970_2030) +
tm_symbols(col = "black", border.col = "white", size = "population_millions") +
tm_facets(by = "year", nrow = 2, free.coords = FALSE)
- another approach would be to create a gif that we can embed
- in this case we use “along” instead of “by”
urb_anim = tm_shape(world) +
tm_polygons() +
tm_shape(urban_agglomerations) +
tm_dots(size = "population_millions") +
tm_facets(along = "year", free.coords = FALSE)
tmap_animation(urb_anim, filename = "urb_anim.gif", delay = 25)interactive maps
- tmap is especially powerful because it allows us to make interactive maps using the same syntax
- all we need to do is enter the interactive mode
tmap_mode("view")## tmap mode set to interactive viewingmap_nzTo go back to regular plotting, we just need enter plotting mode
tmap_mode("plot")## tmap mode set to plotting